Multi-Species Boids For Document Clustering

For some background information on Reynold's Boids algorithm please look at Reynolds Flocking Simulation with Evolution and Genetic Algorithm

Swarm intelligence is used in many different field of studies not only for simulation purposes but also to solve NP-Hard problems sub-optimally. One such example can be found in this beautifully explained video by Sebastian Lague in which he solves the traveling salesman problem using ant swarm intelligence. Another example to showcase the brilliance of swarm intelligence for solving NP-Hard problems is using Reynold's flocking boids simulation to solve document clustering problem.

Document clustering involves the use of descriptors to represent the contents of the document within the cluster. The descriptors are then used to calculate the similarity between the documents. The similarity between the documents is then used to group the documents into clusters. Document clustering is used is many different applications involving gathering of similar data for example getting similar articles relating to a query. The algorithm used to solve the document clustering problem is called the K-Means clustering algorithm. The algorithm takes an iterative approach that starts with a random set of clusters and then iteratively moves the clusters to the center of the documents in the cluster. The algorithm stops when the clusters stop moving. A downside of this algorithm is that it is NP-Hard even though it is complete. This means that for a large enough dataset the algorithm will take an enormous amount of time.

This Paper by Xiaohui Cui, Jesse St. Charles and Thomas Potok explores the idea of utilizing GPU parallelism and boids simulation to solve the document clustering problem. The paper proposes a new algorithm called the Multi-Species Flocking (MSF) algorithm however, since the algorithm is O(n^2), GPU parallelism is used to make the algorithm run a lot faster.

The base code for boids is done in Unity by Sebastian Lague and explained well in this Video. I added the code for the Multi-Species Flocking algorithm by assigning a type to each boid in real time. The type is assigned based on the number of clusters. The boids are then grouped into clusters based on the type of the boid. Different types are represented using different colors. In addition to the MSF I also added CPU level parallelism using unity's job system to get extra performance.

The number of clusters can be changed in real time using the boid manager in the Scene. The algorithm currently lacks target cluster points which causes the boids to move in random directions.

Compute Shader

#pragma kernel CSMain
static const int threadGroupSize = 1024;

struct Boid {
    float3 position;
    float3 direction;

    float3 flockHeading;
    float3 flockCentre;
    float3 separationHeading;
    int numFlockmates;

    int boidType;
};

RWStructuredBuffer(Boid) boids;
int numBoids;
float viewRadius;
float avoidRadius;

[numthreads(threadGroupSize,1,1)]
/* id of the boid we are looking at which loops over all the other boids
 * check the distance between boids and check the type, only update the boids flocks if they are of the same type
 * and the distance is less then the view radius.
 * if the boid we are looking at is of a different type than the boid we are comparing with then add to the seperation force
 * this allows boids of different types to never flock with each other.
 */
void CSMain (uint3 id : SV_DispatchThreadID)
{
    for (int indexB = 0; indexB < numBoids; indexB ++) {
        if (id.x != indexB) {
            Boid boidB = boids[indexB];
            float3 offset = boidB.position - boids[id.x].position;
            float sqrDst = offset.x * offset.x + offset.y * offset.y + offset.z * offset.z;

            if (sqrDst < viewRadius * viewRadius && boids[id.x].boidType == boidB.boidType) {
                boids[id.x].numFlockmates += 1;
                boids[id.x].flockHeading += boidB.direction;
                boids[id.x].flockCentre += boidB.position;

                if (sqrDst < avoidRadius * avoidRadius) {
                    boids[id.x].separationHeading -= offset / sqrDst;
                }
            }
            if(sqrDst < viewRadius * viewRadius && boids[id.x].boidType != boidB.boidType)
                boids[id.x].separationHeading -= offset / sqrDst;
        }
    }
}

Unity Jobs to update boids

[BurstCompile]
public struct BoidUpdateJob : IJobParallelFor
{
    public NativeArray(float3) position;
    public NativeArray(float3) forward;
    public NativeArray(float3) velocity;
    public NativeArray(float3) acceleration;
    public NativeArray(float3) avgFlockHeading;
    public NativeArray(float3) avgAvoidanceHeading;
    public NativeArray(float3) centreOfFlockmates;
    public NativeArray(int) numPerceivedFlockmates;
    [ReadOnly] public float deltaTime;

    public void Execute(int i)
    {
        acceleration[i] = float3.zero;

        if (numPerceivedFlockmates[i] != 0)
        {
            centreOfFlockmates[i] /= numPerceivedFlockmates[i];

            float3 offsetToFlockmatesCentre = (centreOfFlockmates[i] - position[i]);

            var alignmentForce = SteerTowards(avgFlockHeading[i], velocity[i]) * 2f;
            var cohesionForce = SteerTowards(offsetToFlockmatesCentre, velocity[i]) * 1f;
            var seperationForce = SteerTowards(avgAvoidanceHeading[i], velocity[i]) * 2.5f;

            acceleration[i] += (float3) alignmentForce;
            acceleration[i] += (float3) cohesionForce;
            acceleration[i] += (float3) seperationForce;
        }

        velocity[i] += acceleration[i] * deltaTime;
        float x = velocity[i].x;
        float y = velocity[i].y;
        float z = velocity[i].z;
        float speed = Mathf.Sqrt((x * x) + (y * y) + (z * z));
        float3 dir = velocity[i] / speed;
        speed = Mathf.Clamp(speed, 5, 8);
        velocity[i] = dir * speed;

        position[i] += velocity[i] * deltaTime;
        forward[i] = dir;
    }

    float3 SteerTowards(Vector3 vector, Vector3 vel)
    {
        Vector3 v = vector.normalized * 8 - vel;
        return Vector3.ClampMagnitude(v, 8);
    }
}